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Based on a phenomenological model with s± or s-wave pairing symmetry, the mixed-state ef- 
fect on quasiparticle interference in iron-based superconductors is investigated by solving large-scale 
Bogoliubov-de Gennes equations based on the Chebyshev polynomial expansion. Taking into ac- 
count the presence of magnetic field, our result for the s± pairing is in qualitative agreement with 
recent scanning tunneling microscopy experiment while for the s-wave pairing, the result is in appar- 
ent contradiction with experimental observations, thus excluding the s-wave pairing. Furthermore, 
we treat the effect of vortices rigorously instead of approximating the vortices as magnetic impurities, 
thus our results are robust and should be more capable of explaining the experimental data. 

PACS numbers: 74.70.Xa, 74.55.+V, 74.25.Op, 74.25. Ha 



Introduction. — The discovery of iron-based supercon- 
ductors fjj has triggered efforts to elucidate the super- 
conducting (SC) pairing mechanism in these materials. 
One hotly debated issue is the symmetry and struc- 
ture of the SC gap. Theoretically, it was initially sug- 
gested that the pairing may be established via spin fluc- 
tuations, leading to the so-called s± pairing symmetry 
(Ak ~ cosk x + cos k y defined in the 2Fe/cell Brillouin 
zone). In this case the SC order parameter (OP) exhibits 
a sign reversal between the hole pockets (around the T 
point) and electron pockets (around the M point) ||. 
Later, s-wave pairing symmetry without sign reversal was 
also proposed to be a possible candidate which is induced 
by orbital fluctuations due to the electron-phonon inter- 
action Q. Experimentally, the results about the pair- 
ing symmetry remain highly controversial as well. For 
example, in Bao.6Ko.4Fe2 As?, an optimally hole-doped 
iron-based superconductor [4[ , the SC gaps measured by 
the angle resolved photoemission spectroscopy can be ap- 
proximately fitted by | Ak | ~ | cos k x + cos k y | , with al- 
most isotropic gaps on all the Fermi surfaces (FS) 
indicating the possible pairing symmetry to be either 
s± or s-wave. Neutron scattering experiments observed 
a resonance peak at Q = (jr, it) below the SC transi- 
tion temperature as predicted by some theoretical 
works assuming s± symmetry Q, thus at first glance they 
seemed to support each other. However, later theoretical 
works suggested that the experimentally observed reso- 
nance peak can also be reproduced by assuming s-wave 
pairing Q. 

Recently, in order to clarify the pairing symmetry, 
Hanaguri et al. used scanning tunneling microscopy 
(STM) to image the quasiparticle interference (QPI) pat- 
terns in the SC state [9[ . They proposed that the relative 
sign of the SC OP can be determined from the magnetic- 
field dependence of quasiparticle scattering amplitudes 
and claimed that their experimental data were only con- 
sistent with the s± scenario but not the s-wave one. Soon 



after, the experimental results were put into question and 
argued instead to arise from the Bragg scattering but not 
due to the QPI because the observed peaks are too sharp 
fioj . On the other hand, the magnetic field will induce 
vortices into the system and lead to the inhomogeneity of 
the pairing OP in real space, thus affecting the QPI pat- 
terns. Theoretical analyzes performed previously have 
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investigated the mixed-state effect on the QPI 
However, in Ref. [ll[, on ly amplitude suppression of the 
OP near the vortex core was considered, without tak- 
ing into account the phase variation. In Ref. [13], the 
mixed-state effect was treated by using quasiclassical ap- 
proximation and the QPI derived in this way was rather 
broad compared to the experimental observation. Fur- 
thermore, in Ref. [li|, only the effect of the Zeeman 
splitting was considered, which should be negligible with 
respect to the mixed-state effect. Thus there still lacks 
direct theoretical confirmation of the mixed-state effect 
on the QPI. 

Therefore in this work we adopt a phenomenological 
model with s± pairing symmetry to study the influence 
of vortices on the QPI by directly solving large-scale 
Bogoliubov-dc Genncs (BdG) equations in real space 
based on the Chebyshev polynomial expansion. For com- 
parison, the problem is also studied for s-wave pairing. In 
this way the mixed-state effect on the QPI can be rigor- 
ously investigated and the results unambiguously support 
s± pairing symmetry in iron-based superconductors. 

Method. — We begin with an effective two-orbital model 
on a two-dimensional lattice [lil ] , with a phenomenolog- 
ical form for the intraorbital pairing terms. The Hamil- 
tonian can be written as 



h = - Yl [kj^ + ^-Vi 
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Here i and j are the site indices while a, /3 = 1,2 are 
the orbital ones, a represents the spin and fi is the 
chemical potential. Then we consider potential scat- 
tering by nonmagnetic impurities through Vi mp with 



being the locations of the impurities. Ay, 



ij,a/3 



v^Saf, (( Cj - | 3.|.Ci Q ;-|-) — (cjfj^Ciai)) * s the intraorbital spin sin- 
glet bond OP, where Vij is the onsitc [i = j] or ncxt- 
nearest-neighbor [i = j ± (x ± y)] attraction which we 
choose to achieve the s-wave or s± pairing symmetry, 
respectively. In the presence of a magnetic field B per- 
pendicular to the plane, the hopping integral can be 



expressed as t lj a/3 = % iai gexp[i^ J\ A(r) • dr], where 
$0 = hc/2e is the SC flux quantum, and A = (—By, 0, 0) 
is the vector potential in the Landau gauge. The hopping 
integrals are 
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a = j8,« = 3±x(y), 
a = P,i = j±(x + y), 
a = 0,i = j± (x-y), 
a / P,i=j±(x±y), 
otherwise. 

(2) 



The Chebyshev polynomials can be written as <f>k (x) 
cos[fc arccoscc] and statisfy 
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Here r, s, 7 = 1, • • • , i with 7 = 4N x N y and (TVj,) 
being the number of lattice sites along £ (y) direction of 
the 2D lattice, a = (75™°* - E™ m )/{2 - e) (e > is a 
small number), 6 = (£™ aa; + E™ m )/2, w = (w - 6)/a, 
£ 7 = (£? 7 - 6)/a and M = (M - 67) /o. 

If we further define the L-dimcnsional vectors e(o) and 
h(o) as e(o) 7 = <5 70 and /i(o) 7 = <5 70 +27v s (A^ = N x N y 
and o = 1, ■ • ■ , 2iV s ), we can express the self-consistent 
parameters as 



dwf(w)d mm (w) 



doj f (w)e(m) T d(uj)e(m) 



4> (x) = l,<f>x(x) = x,4> k (x) = 2x<fi k - 1 (x) - (f) k - 2 (x), 



f^^-4> k (x)4> k {x') = 5{x-x'), 

, n v k 



(3) 



fe=0 



where W(x) = \j\f\ — x 2 , v k = 7r(l + S k o)/2 and x G 
[—1,1]. Next we define the Green's function matrix: 



V g k e(m) T e k (m) — , 

k— 



2V-1 



n 3 -«4. = 1 - V* g k h(m) T h k (m) — , 
— ffc 



G(r)=-(TTC(r)Ct(0)), 
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with C f = (■ ■ • , ct lf , ct 2f , ■ • ■ , Cjii, c 3 2i, ■■■)■ Eq. (1) 
can be diagonalized as 

H = CHlC = C^QQ^MQQ^C = &D<f>. (5) 

Here Q is a unitary matrix that satisfies (Q'MQ) rs = 
D rs = 8 rs E s and the spectral function can be expressed 
as [H 



d rs (oj) = --—{G rs (uj + irj) - G rs (ui - irj)] 
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where m = 2(j y + N y j x ) + (3 and n = 2(i y + N y i x ) + a 
with i x ,j x = 0, ■ • • , N x - 1 and i y , j y = 0, ■ ■ ■ , N v — 1. 
e k (m) = 4> k (M)e(m), h k (m) = <f> k (M)h(m) and T k = 
f , dujf(auj + b)W(Cj)cj) k {Cj). N is the cutoff in the sum- 
mation and Ofc is the kernel we convolute to avoid Gibbs 
oscillations [l6(. In addition, fix) is the Fermi distribu- 
tion function. At zero temperature we have 



arccos(— -) k = 0, 



n 
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k + 0. 



Then we can solve the BdG equations self-consistently 
and the chemical potential is determined by the doping 
concentration. The calculation is repeated until the ab- 
solute error of the OP between two consecutive iteration 
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steps and that of the total electron number are less than 
10 . The local density of states (LDOS) is given by 

1 N ~ 1 

pM = - J2 E f[W{Cj)<j> k {Q)e{m) T e k {m) 
a p k=0 k 



+W{w')cp k {u')h{m) T h k {m)] 



(9) 



where u>' = (—lu — b)/a. 

The benefits of this method are threefold. First, it 
requires much less storage than the exact diagonaliza- 
tion method since the matrix M is sparse, thus we can 
solve large-scale BdG equations and obtain the QPI by 
Fourier transforming the real-space LDOS in sufficiently 
wide range. Second, it is applicable in parallel computa- 
tion because the self-consistent parameters on each lat- 
tice site can be calculated separately. Third, the expan- 
sion scheme is very stable and efficient. 

In our calculation, the magnitudes of the parame- 
ters are chosen as ti-4 = 1,0.4,-2,0.04. Magnetic 
unit cells are introduced where each unit cell accommo- 
dates four SC flux quanta and the linear dimension is 
N x x N y = 80 x 80, corresponding to a magnetic field 
B w 8.32 Tesla, close to the experimental value (10 Tesla) 
Va and Vij [i = j± (x±y)] are chosen to be —2.8 and 
—2, respectively Moreover we introduce 12 randomly 
distributed impurities with V imp = 0.3. E™ ax (E™ m ) is 
chosen as 1.5 (—1.5) band width. Throughout the paper, 
we set the system to be 20% hole-doped. In calculating 
the self-consistent parameters, we use the Jackson kernel 
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(N - k + 1) cos ^ + sin ^ cot jfc 
N + l 



(10) 



with e = 0.001 and N = 500. For the LDOS we convolute 
the Lorcntz kernel 
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sinh[A(l - |0 
sinh A 
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with A = 4, e = 0.004 and N = A/e. 

Results and discussion. — First we consider the s± case. 



For V, 



imp 



0, Fig- HI a ) shows that there exists 



negative-energy in-gap peak in the LDOS at the vortex 
core center, in agreement with the experimental observa- 
tion [I?} and our previous results based on the exact diag- 
onalization method [ljl], indicating the reliability of the 
Chebyshev polynomial-expansion scheme. Figure HJc) 
plots the spatial distribution of Zj(uj) = pj(u>)/ pj(—u>) 
as defined in Ref. @, for B = and at u = A. 
The locations of the impurities can be clearly identi- 
fied as the low-intensity spots. The corresponding Z q (w), 
which is the Fourier transformation of Zj(oj), is shown in 
Fig. H{d). We notice that there are high-intensity peaks 
at q2 [(±7r, 0) and (0, ±7r)] which arise from the inter- 
pocket scattering between the hole and electron pockets. 
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FIG. 1: (color online) (a) The LDOS at the vortex core center 
(red solid) and at B = (black dot) as a function of the 
reduced energy uu/A. Here 2A is the SC gap between two SC 
coherence peaks in the LDOS at B = 0. The gray-dotted line 
indicates the position of u> — 0. (b) Zq(B 7^ 0) — Z^{B = 0) 
at ui — A. (c) and (d) show Zj(uj = A) and Z^(lo = A) at 
B = 0, respectively, (e) and (f) are the same as (c) and (d), 
respectively, but for B 7^ 0. (c) and (e) share the same color 
scale while the case is similar for (d) and (f). 



Furthermore the intensities at q3 [(±7r, ±7t)] are much 
weaker than those at q2 , consistent with the experimental 
observation Upon applying the magnetic field, vor- 
tices are introduced into the system and their locations 
are denoted as the high- intensity spots in Fig. [TJe). In 
this case, from Fig. [TJf ) we can see that there exist addi- 
tional peaks at q3 , whose intensities are enhanced by the 
application of the magnetic field and they are due to the 
interpocket scattering between different electron pockets 
(see Fig. 1 in Ref. 0). Figure [1Tb) shows the mag- 
netic field-induced change in the QPI intensities defined 
as Zq(B 7^ 0) — Zq(B = 0). In the presence of the time- 
reversal symmetry breaking due to the magnetic field, 
the sign-preserving scatterings at q3 (q3 connects the FS 
with the same sign of the SC OP) are enhanced while the 
sign-reversing scatterings at q2 (q2 connects the FS with 
the opposite sign of the SC OP) are suppressed. Both 
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FIG. 2: (color online) (a), (b), (c) and (d) are similar to Figs. 
Ha), [T]L), QJd) and[TJf), respectively, (c) and (d) share the 
same color scale. 

the locations and the sharpness of the QPI peaks shown 
in Figs. (Ub),[IJd) andQJf) are in reasonable agreement 
with the experimental data suggesting that the ex- 
perimentally observed peaks are indeed due to the QPI 
but not the Bragg scattering as argued in Ref. . 

Next we consider the s-wave case. For Vi mp = 0, the 
LDOS at B = (black dot) and at the vortex core center 
(red solid) shown in Fig. [2ja) are also consistent with 
those obtained by exact diagonalization [3], again sug- 
gesting the validity of the current polynomial-expansion 
scheme. From Fig. [2^c) we notice, the QPI in the pres- 
ence of impurities at B = exhibits no pronounced peaks 
at either q 2 or q 3 , as compared to the clear peaks at q 2 
in the s± case as shown in Fig. [TJd). After applying the 
magnetic field, the intensities at q3 are enhanced and 
they form sharp peaks as shown in Fig. GJd), similar 
to the s± case. At last, from the magnetic field-induced 
change in the QPI intensities plotted in Fig. H^b) we can 
see, the intensities are enhanced at q 3 but remain almost 
unchanged at q 2 . The lack of distinct structures at q 2 
is in stark contrast to the s± case and is inconsistent 
with the experimental observations ||. Therefore, the 
different behavior of the QPI intensities at q 2 in the s± 
and s-wave pairing cases makes it possible to distinguish 
these two types of pairing symmetry since the STM ex- 
periment observed clear structures at q 2 , thus excluding 
the possibility of s-wave pairing in iron-based supercon- 
ductors. 

Summary. — In summary, by using the Chebyshev poly- 
nomial expansion to directly solve large-scale BdG equa- 



tions in real space, we have investigated the mixed-state 
effect on QPI in iron-based superconductors by assum- 
ing s± or s-wave pairing symmetry. For the s± pairing, 
the QPI intensities at q 2 which connects the FS with 
the opposite sign of the SC OP are suppressed by the 
application of the magnetic field while the situation at 
q 3 is reversed where q3 connects the FS with the same 
sign of the SC OP. The obtained results at both B = 
and B ^ are in qualitative agreement with experiment, 
suggesting that the experimentally observed peaks are in- 
deed due to QPI. On the other hand, for the s-wave pair- 
ing, the QPI intensities at q 2 arc featureless both with 
and without the magnetic field. Based on the available 
experimental data, the s-wave pairing can be excluded in 
iron-based superconductors. 
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